
###########################################################
##### Haiti elite network project  		          			#####
##### fam level regressions	                      		#####
##### 2021 mar 03                   									#####
###########################################################


#####
## read in data
#####

fam <- read.csv('01_Data/02_Clean/fam.csv')
all <- read.csv('01_Data/02_Clean/allfams.csv')


#####
## summary stats
#####

##
## check coverage of gen data
##

prop.table(table(is.na(all$bonw_02_wnind[all$pol==1])))
prop.table(table(is.na(all$bonw_02_wnind[all$mil==1])))
prop.table(table(is.na(all$bonw_02_wnind[all$biz==1])))


##
## TABLE A2 - Summary stats importers only
##

fam <- data.table(fam)
fam$value_1_mil <- fam$value_1/1000000

temp1 <- fam[,lapply(.SD, FUN = function(x) mean(x, na.rm = T)), by = 'coup',
             .SDcols = c('immig', 'syrian', 'pol', 'mil', 'bonw_02_wnind', 'degree_all_uw',
                         'nind', 'reachability',  'cshare', 'noconsump',
                         'wgt_pct', 'value_1_mil',
                         'bulk_ln', 'divis_ln', 'ref_lib', 'time_con', 'pci_value')]
temp2 <- fam[,lapply(.SD, FUN = function(x) sd(x, na.rm = T)), by = 'coup',
             .SDcols = c('immig', 'syrian', 'pol', 'mil', 'bonw_02_wnind', 'degree_all_uw',
                         'nind', 'reachability',
                         'wgt_pct', 'value_1_mil',  'cshare', 'noconsump',
                         'bulk_ln', 'divis_ln', 'ref_lib', 'time_con', 'pci_value')]
temp3 <- fam[,lapply(.SD, FUN = function(x) length(na.omit(x))), by = 'coup',
             .SDcols = c('immig', 'syrian', 'pol', 'mil', 'bonw_02_wnind', 'degree_all_uw',
                         'nind', 'reachability',
                         'wgt_pct', 'value_1_mil', 'cshare', 'noconsump',
                         'bulk_ln', 'divis_ln', 'ref_lib', 'time_con', 'pci_value')]

tab <- rbind(temp3[1,], temp1[1,], temp2[1,], temp3[2,], temp1[2,], temp2[2,])
tab <- t(tab)
rownames(tab) <- c('Coup', 'Immigrant', 'Middle Eastern', 'Political elite', 'Military elite',
                   'Bonacich centrality', 'Degree', 'Family size', 'Reachability',
                   'Market share', 'Value (mil USD)', 'Consumption share', 
                   'All inputs', 'Bulkiness', 'Divisibility', 'Reference price',
                   'Time sensitivity', 'Complexity')
xtable(tab, digits = c(0, 0, 2, 2, 0, 2, 2))


##
## TABLE A3 - Summary stats all elite
##

all <- data.table(all)
temp4 <- all[,lapply(.SD, FUN = function(x) mean(x, na.rm = T)), by = 'coup',
             .SDcols = c('immig', 'syrian', 'pol', 'mil', 'biz', 'bonw_02_wnind', 'degree_all_uw',
                         'nind', 'reachability')]
temp5 <- all[,lapply(.SD, FUN = function(x) sd(x, na.rm = T)), by = 'coup',
             .SDcols = c('immig', 'syrian', 'pol', 'mil', 'biz', 'bonw_02_wnind', 'degree_all_uw',
                         'nind', 'reachability')]
temp6 <- all[,lapply(.SD, FUN = function(x) length(na.omit(x))), by = 'coup', 
             .SDcols = c('immig', 'syrian', 'pol', 'mil', 'biz', 'bonw_02_wnind', 'degree_all_uw',
                         'nind', 'reachability')]

tab <- rbind(temp6[2,], temp4[2,], temp5[2,], temp6[1,], temp4[1,], temp5[1,])
tab <- t(tab)
rownames(tab) <- c('Coup', 'Immigrant', 'Middle Eastern', 'Political elite', 'Military elite',
                   'Business elite',
                   'Bonacich centrality', 'Degree', 'Family size', 'Reachability')
xtable(tab, digits = c(0, 0, 2, 2, 0, 2, 2))



#####
## Figure 2: plot coup participation for deciles of centrality
#####

## deciles

## importers sample 

fam$bonw_02_wnind_dec <- cut(fam$bonw_02_wnind_st, 
                             breaks = quantile(fam$bonw_02_wnind_st, seq(0,1,0.1), na.rm=T),
                             labels = paste0('Q', seq(1,10,1)))

t <- melt(fam, id.vars = c('bonw_02_wnind_dec'), measure.vars = c('coup'))
t1 <- cast(t, bonw_02_wnind_dec ~ variable, mean)
t2 <- cast(t, bonw_02_wnind_dec ~ variable, function(x) sd(x)/sqrt(length(x)))
setnames(t1, 'coup', 'mean_coup')
setnames(t2, 'coup', 'se_coup')
t <- merge(t1, t2, by = 'bonw_02_wnind_dec')
t <- subset(t, is.na(t$bonw_02_wnind_dec)==F)

pdf('03_Figures/fam_partic_bin.pdf')
ggplot(t, aes(x = bonw_02_wnind_dec, y = mean_coup, group = 1)) +
  geom_line() +
  geom_point() +
  geom_ribbon(aes(ymin = mean_coup - 1.96*se_coup, 
                  ymax = mean_coup + 1.96*se_coup), alpha = 0.1) +
  ylim(0,1) +
  theme_minimal() +
  theme(text = element_text(size=16), 
        axis.text.x = element_text(angle=45, vjust=0.95, hjust=1),
        legend.position = "none",
        panel.spacing = unit(0.5, "lines")) +
  labs(y = "Coup Participation", 
       x = "Centrality")
dev.off()


## all elite sample

all$bonw_02_wnind_dec <- cut(all$bonw_02_wnind_st, 
                             breaks = quantile(all$bonw_02_wnind_st, seq(0,1,0.1), na.rm=T),
                             labels = paste0('Q', seq(1,10,1)))

t <- melt(all, id.vars = c('bonw_02_wnind_dec'), measure.vars = c('coup'))
t1 <- cast(t, bonw_02_wnind_dec ~ variable, mean)
t2 <- cast(t, bonw_02_wnind_dec ~ variable, function(x) sd(x)/sqrt(length(x)))
setnames(t1, 'coup', 'mean_coup')
setnames(t2, 'coup', 'se_coup')
t <- merge(t1, t2, by = 'bonw_02_wnind_dec')
t <- subset(t, is.na(t$bonw_02_wnind_dec)==F)

pdf('03_Figures/all_partic_bin.pdf')
ggplot(t, aes(x = bonw_02_wnind_dec, y = mean_coup, group = 1)) +
  geom_line() +
  geom_point() +
  geom_ribbon(aes(ymin = mean_coup - 1.96*se_coup, 
                  ymax = mean_coup + 1.96*se_coup), alpha = 0.1) +
  ylim(0,1) +
  theme_minimal() +
  theme(text = element_text(size=16), 
        axis.text.x = element_text(angle=45, vjust=0.95, hjust=1),
        legend.position = "none",
        panel.spacing = unit(0.5, "lines")) +
  labs(y = "Coup Participation", 
       x = "Centrality")
dev.off()




##
## TABLE 1: Coup participation
##

all$clust <-  as.factor(all$fastgreedy)
fam$clust <-  as.factor(fam$fastgreedy)


mod4 <- lm(coup ~ bonw_02_wnind_st, data = fam)
mod4$se <- coeftest(mod4, vcov = vcovHC, type = "HC1")[,2]
mod4$se_clust <- coeftest(mod4, vcov = cluster.vcov(mod4, fam$clust))[,2]
inc.obs <- complete.cases(fam[,colnames(mod4$model)])
mod4$nclust <- length(unique(fam$clust[inc.obs==T]))

mod5 <- lm(coup ~ bonw_02_wnind_st + nind_log, data = fam)
mod5$se <- coeftest(mod5, vcov = vcovHC, type = "HC1")[,2]
mod5$se_clust <- coeftest(mod5, vcov = cluster.vcov(mod5, fam$clust))[,2]
inc.obs <- complete.cases(fam[,colnames(mod5$model)])
mod5$nclust <- length(unique(fam$clust[inc.obs==T]))

mod6 <- lm(coup ~ bonw_02_wnind_st + value_log_st + cshare2_st + noconsump + nind_log,
           data = fam)
mod6$se <- coeftest(mod6, vcov = vcovHC, type = "HC1")[,2]
mod6$se_clust <- coeftest(mod6, vcov = cluster.vcov(mod6, fam$clust))[,2]
inc.obs <- complete.cases(fam[,colnames(mod6$model)])
mod6$nclust <- length(unique(fam$clust[inc.obs==T]))

mod7 <- lm(coup ~ bonw_02_wnind_st + value_log_st + cshare2_st + noconsump +
             pol + mil + syrian + immig + nind_log,
           data = fam)
mod7$se <- coeftest(mod7, vcov = vcovHC, type = "HC1")[,2]
mod7$se_clust <- coeftest(mod7, vcov = cluster.vcov(mod7, fam$clust))[,2]
inc.obs <- complete.cases(fam[,colnames(mod7$model)])
mod7$nclust <- length(unique(fam$clust[inc.obs==T]))

mod8 <- lm(coup ~ bonw_02_wnind_st + value_log_st + cshare2_st + noconsump +
             ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st + 
             pol + mil + syrian + immig + nind_log, data = fam)
mod8$se <- coeftest(mod8, vcov = vcovHC, type = "HC1")[,2]
mod8$se_clust <- coeftest(mod8, vcov = cluster.vcov(mod8, fam$clust))[,2]
inc.obs <- complete.cases(fam[,colnames(mod8$model)])
mod8$nclust <- length(unique(fam$clust[inc.obs==T]))

mod8b <- lm(coup ~ bonw_02_wnind_st + value_log_st + cshare2_st + noconsump +
             ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st + 
             pol + mil + syrian + immig + nind_log + clust, data = fam)
mod8b$se <- coeftest(mod8b, vcov = vcovHC, type = "HC1")[,2]
mod8b$se_clust <- coeftest(mod8b, vcov = cluster.vcov(mod8b, fam$clust))[,2]
inc.obs <- complete.cases(fam[,colnames(mod8b$model)])
mod8b$nclust <- length(unique(fam$clust[inc.obs==T]))

mod1 <- lm(coup ~ bonw_02_wnind_st, data = all)
mod1$se <- coeftest(mod1, vcov = vcovHC, type = "HC1")[,2]
mod1$se_clust <- coeftest(mod1, vcov = cluster.vcov(mod1, all$clust))[,2]
inc.obs <- complete.cases(all[,colnames(mod1$model)])
mod1$nclust <- length(unique(all$clust[inc.obs==T]))

mod2 <- lm(coup ~ bonw_02_wnind_st + nind_log, data = all)
mod2$se <- coeftest(mod2, vcov = vcovHC, type = "HC1")[,2]
mod2$se_clust <- coeftest(mod2, vcov = cluster.vcov(mod2, all$clust))[,2]
inc.obs <- complete.cases(all[,colnames(mod2$model)])
mod2$nclust <- length(unique(all$clust[inc.obs==T]))

mod3 <- lm(coup ~ bonw_02_wnind_st + biz + pol + mil + syrian + immig + nind_log, data = all)
mod3$se <- coeftest(mod3, vcov = vcovHC, type = "HC1")[,2]
mod3$se_clust <- coeftest(mod3, vcov = cluster.vcov(mod3, all$clust))[,2]
inc.obs <- complete.cases(all[,colnames(mod3$model)])
mod3$nclust <- length(unique(all$clust[inc.obs==T]))

mod3b <- lm(coup ~ bonw_02_wnind_st + biz + pol + mil + syrian + immig + nind_log + clust, data = all)
mod3b$se <- coeftest(mod3b, vcov = vcovHC, type = "HC1")[,2]
mod3b$se_clust <- coeftest(mod3b, vcov = cluster.vcov(mod3b, all$clust))[,2]
inc.obs <- complete.cases(all[,colnames(mod3b$model)])
mod3b$nclust <- length(unique(all$clust[inc.obs==T]))


## for paper

stargazer(mod4, mod5, mod6, mod7, mod8, mod8b, mod1, mod2, mod3, mod3b, 
          covariate.labels = c("Centrality", "Family Size",
                               "Political Elite", "Military Elite", "Business Elite",
                               "Middle Eastern", "Immigrant",
                               "Value", 'Consumption Share', 'All Inputs',
                               "Reference Price", "Complexity", "Divisibility",
                               "Bulkiness", "Time Sensitivity",
                               "Constant"),
          se = list(mod4$se, 
                    mod5$se, mod6$se, mod7$se, mod8$se, mod8b$se,
                    mod1$se, mod2$se, mod3$se, mod3b$se),
         order = c('bonw', 'nind_log', 'pol', 'mil', 'biz', 'syrian', 'immig', 'value_log', 'cshare', 'nocons',
                   'ref_lib', 'pci_value', 'divis_ln', 'bulk_ln', 'time_con'),
           omit = c('clust'),
          add.lines = list(c('Community FE', '', '', '', '\\checkmark', '', '', '', '', '', '\\checkmark'),
                           c('Product Characteristics', '', '', '', '', '', '', '', '', '\\checkmark', '\\checkmark'),
                           c('Clusters', round(c(mod4$nclust, mod5$nclust, mod6$nclust, mod7$nclust, mod8$nclust, mod8b$nclust, mod1$nclust, mod2$nclust, mod3$nclust, mod3b$nclust), digits = 1)),
                           c('Clustered SEs', round(c(mod4$se_clust, mod5$se_clust, mod6$se_clust, mod7$se_clust, mod8$se_clust, mod8b$se_clust, mod1$se_clust, mod2$se_clust, mod3$se_clust, mod3b$se_clust), digits = 3)),
                           c('Sample', '\\multicolumn{6}{c}{Importers}', '\\multicolumn{4}{c}{All elite}')),
          no.space = T, keep.stat = c('n', 'rsq'), digits = 3)


## predicted probs

predict(mod8b, newdata = data.frame('bonw_02_wnind_st' = c(0,1), 'biz' = rep(mean(fam$biz,na.rm=T),2), 
                                   'pol' = rep(mean(fam$pol,na.rm=T),2), 'mil' = rep(mean(fam$mil,na.rm=T),2),
                                   'syrian' = rep(mean(fam$syrian,na.rm=T),2), 'immig' = rep(mean(fam$biz,na.rm=T),2),
                                   'noconsump' = rep(mean(fam$noconsump,na.rm=T),2), 
                                   'ref_lib_st' = rep(mean(fam$ref_lib_st,na.rm=T),2),
                                   'pci_value_st' = rep(mean(fam$pci_value_st,na.rm=T),2),
                                   'divis_ln_st' = rep(mean(fam$divis_ln_st,na.rm=T),2),
                                   'bulk_ln_st' = rep(mean(fam$bulk_ln_st,na.rm=T),2),
                                   'time_con_st' = rep(mean(fam$time_con_st,na.rm=T),2),
                                   'value_log_st' = rep(mean(fam$value_log_st,na.rm=T),2),
                                   'cshare2_st' = rep(mean(fam$cshare2_st,na.rm=T),2),
                                   'nind_log'=rep(mean(fam$nind_log,na.rm=T),2),
                                   'clust'=rep(get_mode(na.omit(fam$clust)),2)))




#####
## Figure 3: robustness to bonacich measures with varying 1/delta
#####

fam <- data.frame(fam)
all <- data.frame(all)

## importers data

coef <- paste0('bonw',seq(1,9,1),'_wnind')
for (i in 1:length(coef)){
  t <- fam[,coef[i]]
  t_st <- (t-mean(t, na.rm=T))/sd(t, na.rm=T)
  fam <- cbind.data.frame(fam, t_st)
  setnames(fam, c('t_st'), c(paste0(coef[i], '_st')))    
}

vars <- c(paste0(rev(coef), '_st'))
betas <- matrix(NA, length(vars), 6)

for (i in 1:length(vars)){
  fam$t <- fam[,vars[i]]
  mod <- lm(coup ~ t, data = fam)
  betas[i,1] <- summary(mod)$coef['t',1]
  betas[i,2] <- coeftest(mod, vcov = vcovHC, type = "HC1")['t',2]
}

betas[,3] <- betas[,1] + 1.64*betas[,2]
betas[,4] <- betas[,1] - 1.64*betas[,2]
betas[,5] <- betas[,1] + 1.96*betas[,2]
betas[,6] <- betas[,1] - 1.96*betas[,2]

pdf('03_Figures/coef_centrangeweight.pdf')
par(xpd = F, mfrow = c(1,1), mar = c(6, 5, 4, 4))
n <- length(vars)
plot(x = NULL, y = NULL, xlim = c(0, dim(betas)[1]+1), ylim = c(-.05, .35), xaxt = "n", 
     ylab = 'Coefficient',
     xlab = "", 
     cex.axis = 1.5, cex.lab = 1.5)
points(x = seq(1, dim(betas)[1], 1), y = betas[,1], col = 'firebrick1', pch = 16, cex = 1.3)
for (i in 1:n){
  points(x = rep(i, 2), y = c(betas[i,3], betas[i,4]), 
         col = 'firebrick1', type = "l", lwd = 4)
  points(x = rep(i, 2), y = c(betas[i,5], betas[i,6]), 
         col = 'firebrick1', type = "l", lwd = 2)
}
abline(h = 0, col = 'black', lty = 2)
mtext(expression(paste("Scaling parameter ", frac(1,delta))),
      cex = 1.2, padj = 14)
axis(side = 1, at = seq(1,n,1), padj = .4,
     # rep("", n),
     labels = c(expression(frac(9,10 * lambda)),
       expression(frac(8,10 * lambda)),
       expression(frac(7,10 * lambda)),
       expression(frac(6,10 * lambda)),
       expression(frac(5,10 * lambda)),
       expression(frac(4,10 * lambda)),
       expression(frac(3,10 * lambda)),
       expression(frac(2,10 * lambda)),
       expression(frac(1,10 * lambda))))
dev.off()


## all elite data

coef <- c(paste0('bonw', seq(1, 9, 1), '_wnind'))
for (i in 1:length(coef)){
  t <- all[,coef[i]]
  t_st <- (t-mean(t, na.rm=T))/sd(t, na.rm=T)
  all <- cbind.data.frame(all, t_st)
  setnames(all, c('t_st'), c(paste0(coef[i], '_st')))    
}

vars <- c(paste0(rev(coef), '_st'))
betas <- matrix(NA, length(vars), 6)

for (i in 1:length(vars)){
  all$t <- all[,vars[i]]
  mod <- lm(coup ~ t, data = all)
  betas[i,1] <- summary(mod)$coef['t',1]
  betas[i,2] <- coeftest(mod, vcov = vcovHC, type = "HC1")['t',2]
}

betas[,3] <- betas[,1] + 1.64*betas[,2]
betas[,4] <- betas[,1] - 1.64*betas[,2]
betas[,5] <- betas[,1] + 1.96*betas[,2]
betas[,6] <- betas[,1] - 1.96*betas[,2]

pdf('03_Figures/coef_centrangeweight_all.pdf')
par(xpd = F, mar = c(5, 5, 4, 4))
n <- length(vars)
plot(x = NULL, y = NULL, xlim = c(0, n+1), ylim = c(-.05, .35), xaxt = "n", ylab = 'Coefficient',
     xlab = "", cex.axis = 1.5, cex.lab = 1.5)
points(x = seq(1, n, 1), y = betas[,1], col = 'firebrick1', pch = 16, cex = 1.3)
for (i in 1:n){
  points(x = rep(i, 2), y = c(betas[i,3], betas[i,4]), 
         col = 'firebrick1', type = "l", lwd = 4)
  points(x = rep(i, 2), y = c(betas[i,5], betas[i,6]), 
         col = 'firebrick1', type = "l", lwd = 2)
}
abline(h = 0, col = 'black', lty = 2)
mtext(expression(paste("Scaling parameter ", frac(1,delta))),
      cex = 1.2, padj = 14)
axis(side = 1, at = seq(1,n,1), padj = .4,
     labels = c(expression(frac(9,10 * lambda)),
                expression(frac(8,10 * lambda)),
                expression(frac(7,10 * lambda)),
                expression(frac(6,10 * lambda)),
                expression(frac(5,10 * lambda)),
                expression(frac(4,10 * lambda)),
                expression(frac(3,10 * lambda)),
                expression(frac(2,10 * lambda)),
                expression(frac(1,10 * lambda))))
dev.off()



#####
## Table A5: Determinants of centrality
#####

mod1 <- lm(bonw_02_wnind_st ~ syrian + immig + reachability + mil + pol + biz, all)
mod1$se <- coeftest(mod1, vcov = vcovHC, type = "HC1")[,2]

mod2 <- lm(bonw_02_wnind_st ~ syrian + immig + reachability + mil + pol + biz + nind_log, all)
mod2$se <- coeftest(mod2, vcov = vcovHC, type = "HC1")[,2]

mod3 <- lm(bonw_02_wnind_st ~ syrian + immig + reachability + mil + pol, fam)
mod3$se <- coeftest(mod3, vcov = vcovHC, type = "HC1")[,2]

mod4 <- lm(bonw_02_wnind_st ~ syrian + immig + reachability + mil + pol + nind_log, fam)
mod4$se <- coeftest(mod4, vcov = vcovHC, type = "HC1")[,2]

mod5 <- lm(bonw_02_wnind_st ~ syrian + immig + reachability + mil + pol + nind_log + value_1_mil, fam)
mod5$se <- coeftest(mod5, vcov = vcovHC, type = "HC1")[,2]

mod6 <- lm(bonw_02_wnind_st ~ syrian + immig + nind_log + reachability + mil + pol + value_1_mil +
             cshare2_st + noconsump + ref_lib_st + pci_value_st + time_con_st + bulk_ln_st + divis_ln_st, fam)
mod6$se <- coeftest(mod6, vcov = vcovHC, type = "HC1")[,2]


stargazer(mod3, mod4, mod5, mod6, mod1, mod2, 
          covariate.labels = c("Middle Eastern", "Immigrant", "Reachability",
                               "Military", 'Political', 'Business',
                               "Family Size (Log)", "Business Value (Mil USD)",
                               "Consumption Share", "All Inputs",
                               "Reference Price", "Complexity",
                               "Time Sensitivity", "Bulkiness", "Divisibility",
                               "Constant"),
          order = c('syrian', 'immig', 'reacha', '^mil$', 'pol', 'biz', 'nind_',
                    'value_1', 'cshare', 'noconsum', 'ref_lib', 'pci_valu',
                    'time_con', 'bulk_ln', 'divis_ln', 'constant'),
          add.lines = list(c('Sample', '\\multicolumn{3}{c}{Importers}', '\\multicolumn{3}{c}{All elite}')),
          no.space = T, keep.stat = c('n', 'rsq'), digits = 2,
          se = list(mod3$se, mod4$se, mod5$se, mod6$se, mod1$se, mod2$se))



#####
## Table A6: Centrality in the Network of Coup Participators 
#####

mod1 <- lm(coup ~ degree_coupall_wnind, data = all)
mod1$se <- coeftest(mod1, vcov = vcovHC, type = "HC1")[,2]

mod2 <- lm(coup ~ degree_coupall_wnind, data = fam)
mod2$se <- coeftest(mod2, vcov = vcovHC, type = "HC1")[,2]


stargazer(mod2, mod1, 
          keep = c('degree_coupall_wnind'),
          covariate.labels = c("Coup Degree"),
          add.lines = list(c('Sample', '\\multicolumn{1}{c}{Importers}', '\\multicolumn{1}{c}{All elite}')),
          no.space = T, keep.stat = c('n', 'rsq'), digits = 3,
          se = list(mod2$se, mod1$se))



#####
## Table A7: weights based on data quality
#####


## make table w weighted regression

all$w <- all$reachability
fam$w <- fam$reachability

all$clust <-  as.factor(all$fastgreedy)
fam$clust <-  as.factor(fam$fastgreedy)

t <- subset(fam, select =c('bonw_02_wnind_st', 'w', 'coup'))
mod4w <- lm(coup ~ bonw_02_wnind_st, data = fam, weights = w)
mod4w$se <- coeftest(mod4w, vcov = vcovHC, type = "HC1")[,2]

mod5w <- lm(coup ~ bonw_02_wnind_st + nind_log, data = fam, weights = w)
mod5w$se <- coeftest(mod5w, vcov = vcovHC, type = "HC1")[,2]

mod6w <- lm(coup ~ bonw_02_wnind_st + value_log_st + cshare2_st + noconsump + nind_log,
            data = fam, weights = w)
mod6w$se <- coeftest(mod6w, vcov = vcovHC, type = "HC1")[,2]

mod7w <- lm(coup ~ bonw_02_wnind_st + value_log_st + cshare2_st + noconsump +
              pol + mil + syrian + immig + nind_log,
            data = fam, weights = w)
mod7w$se <- coeftest(mod7w, vcov = vcovHC, type = "HC1")[,2]

mod8w <- lm(coup ~ bonw_02_wnind_st + value_log_st + cshare2_st + noconsump +
              ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st +
              pol + mil + syrian + immig + nind_log, data = fam, weights = w)
mod8w$se <- coeftest(mod8w, vcov = vcovHC, type = "HC1")[,2]

mod8bw <- lm(coup ~ bonw_02_wnind_st + value_log_st + cshare2_st + noconsump +
               ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st +
               pol + mil + syrian + immig + nind_log + clust, data = fam, weights = w)
mod8bw$se <- coeftest(mod8bw, vcov = vcovHC, type = "HC1")[,2]
inc.obs <- complete.cases(fam[,c(colnames(mod8bw$model)[1:16],'w')])
mod8bw$nclust <- length(unique(fam$clust[inc.obs==T]))

mod1w <- lm(coup ~ bonw_02_wnind_st, data = all, weights = w)
mod1w$se <- coeftest(mod1w, vcov = vcovHC, type = "HC1")[,2]

mod2w <- lm(coup ~ bonw_02_wnind_st + nind_log, data = all, weights = w)
mod2w$se <- coeftest(mod2w, vcov = vcovHC, type = "HC1")[,2]

mod3w <- lm(coup ~ bonw_02_wnind_st + biz + pol + mil + syrian + immig + nind_log + clust,
            data = all, weights = w)
mod3w$se <- coeftest(mod3w, vcov = vcovHC, type = "HC1")[,2]

mod3bw <- lm(coup ~ bonw_02_wnind_st + biz + pol + mil + syrian + immig + nind_log + clust,
             data = all, weights = w)
mod3bw$se <- coeftest(mod3bw, vcov = vcovHC, type = "HC1")[,2]
inc.obs <- complete.cases(fam[,c(colnames(mod3bw$model)[1:9],'w')])
mod3bw$nclust <- length(unique(all$clust[inc.obs==T]))


stargazer(mod4w, mod5w, mod6w, mod7w, mod8w, mod8bw, mod1w, mod2w, mod3w, mod3bw, 
          covariate.labels = c("Centrality", "Family Size"),
          keep = c('bonw_', 'nind_log'),
          omit = c('clust'),
          add.lines = list(c('Community FE', '', '', '', '', '', '\\checkmark', '', '', '', '\\checkmark'),
                           c('Social Characteristics', '', '\\checkmark', '\\checkmark', '\\checkmark', '', '', '', '\\checkmark', '\\checkmark', '\\checkmark'),
                           c('Economic Characteristics', '', '\\checkmark', '\\checkmark', '\\checkmark', '\\checkmark', '', '', '', '', ''),
                           c('Product Characteristics', '', '', '', '', '\\checkmark', '\\checkmark', '', '', '', ''),
                           c('Community FE', '', '', '', '', '', round(mod8bw$nclust,1), '', '', '', round(mod3bw$nclust,1)),
                           c('Sample', '\\multicolumn{6}{c}{Importers}', '\\multicolumn{4}{c}{All elite}')),
          se = list(mod4w$se,
                    mod5w$se, mod6w$se, mod7w$se, mod8w$se, mod8bw$se,
                    mod1w$se, mod2w$se, mod3w$se, mod3bw$se),
          no.space = T, keep.stat = c('n', 'rsq'), digits = 3)



#####
## TABLE A8: Robustness to varying node and edge weights in the centrality measure
#####

fam$bonw_02_wnind_own0_st <- standardize(fam$bonw_02_wnind_own0)
all$bonw_02_wnind_own0_st <- standardize(all$bonw_02_wnind_own0)

vars <- c('bonw_02_wnind_st', 'bonw_02_uw_st',
          'bon_02_wnind_st', 'bon_02_uw_st',
          'bonw0212_02_wnind_st', 'bonwautoc_02_wnind_st', 
          'bonwbin_02_wnind_st', 'bonw_02_wnind2_st',
          'bonwp_02_wnind_st',
          'bonwb_02_wnind_st', 'bonwr_02_wnind_st',
          'bonw_02_wnind_own0_st')

all <- as.data.frame(all); fam <- as.data.frame(fam)

coefs <- matrix(NA, 15, 10)
ses <- matrix(NA, 15, 10)
ps <- matrix(NA, 15, 10)

for (i in 1:length(vars)){
  all$bon_st <- all[,vars[i]]
  fam$bon_st <- fam[,vars[i]]
  # run models
  mod1 <- lm(coup ~ bon_st, all)
  mod2 <- lm(coup ~ bon_st + nind_log, all)
  mod3 <- lm(coup ~ bon_st + biz + pol + mil + syrian + immig + nind_log, all)
  mod3b <- lm(coup ~ bon_st + biz + pol + mil + syrian + immig + nind_log + clust, all)
  mod4 <- lm(coup ~ bon_st, fam)
  mod5 <- lm(coup ~ bon_st + value_log_st + cshare2_st + noconsump, fam)
  mod6 <- lm(coup ~ bon_st + value_log_st + cshare2_st + noconsump +
               ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st, fam)
  mod7 <- lm(coup ~ bon_st + value_log_st + cshare2_st + noconsump +
               ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st +
               biz + pol + mil + syrian + immig, fam)
  mod8 <- lm(coup ~ bon_st + value_log_st + cshare2_st + noconsump +
               ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st +
               pol + mil + syrian + immig + nind_log, fam)
  mod8b <- lm(coup ~ bon_st + value_log_st + cshare2_st + noconsump +
                ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st +
                pol + mil + syrian + immig + nind_log + clust, fam)
  mods <- list(mod4, mod5, mod6, mod7, mod8, mod8b, mod1, mod2, mod3, mod3b)
  for (t in 1:length(mods)){
    mod <- mods[[t]]
    se <- coeftest(mod, vcov = vcovHC, type = "HC1")
    coefs[i,t] <- summary(mod)$coef['bon_st',1]
    ses[i,t] <- se['bon_st',2]
    ps[i,t] <- se['bon_st',4]
  }
}

digits = 2
tab <- rbind(paste0(round(coefs[1,],digits), 
                    ifelse(ps[1,]<0.01, "***", ifelse(ps[1,]<0.05, "**", ifelse(ps[1,]<0.1, "*", "")))),
             paste0("(",round(ses[1,],digits),")"),
             paste0(round(coefs[2,],digits), 
                    ifelse(ps[2,]<0.01, "***", ifelse(ps[2,]<0.05, "**", ifelse(ps[2,]<0.1, "*", "")))),
             paste0("(",round(ses[2,],digits),")"),
             paste0(round(coefs[3,],digits), 
                    ifelse(ps[3,]<0.01, "***", ifelse(ps[3,]<0.05, "**", ifelse(ps[3,]<0.1, "*", "")))),
             paste0("(",round(ses[3,],digits),")"),
             paste0(round(coefs[4,],digits), 
                    ifelse(ps[4,]<0.01, "***", ifelse(ps[4,]<0.05, "**", ifelse(ps[4,]<0.1, "*", "")))),
             paste0("(",round(ses[4,],digits),")"),
             paste0(round(coefs[5,],digits), 
                    ifelse(ps[5,]<0.01, "***", ifelse(ps[5,]<0.05, "**", ifelse(ps[5,]<0.1, "*", "")))),
             paste0("(",round(ses[5,],digits),")"),
             paste0(round(coefs[6,],digits), 
                    ifelse(ps[6,]<0.01, "***", ifelse(ps[6,]<0.05, "**", ifelse(ps[6,]<0.1, "*", "")))),
             paste0("(",round(ses[6,],digits),")"),
             paste0(round(coefs[7,],digits), 
                    ifelse(ps[7,]<0.01, "***", ifelse(ps[7,]<0.05, "**", ifelse(ps[7,]<0.1, "*", "")))),
             paste0("(",round(ses[7,],digits),")"),
             paste0(round(coefs[8,],digits), 
                    ifelse(ps[8,]<0.01, "***", ifelse(ps[8,]<0.05, "**", ifelse(ps[8,]<0.1, "*", "")))),
             paste0("(",round(ses[8,],digits),")"),
             paste0(round(coefs[9,],digits), 
                    ifelse(ps[9,]<0.01, "***", ifelse(ps[9,]<0.05, "**", ifelse(ps[9,]<0.1, "*", "")))),
             paste0("(",round(ses[9,],digits),")"),
             paste0(round(coefs[10,],digits), 
                    ifelse(ps[10,]<0.01, "***", ifelse(ps[10,]<0.05, "**", ifelse(ps[10,]<0.1, "*", "")))),
             paste0("(",round(ses[10,],digits),")"),
             paste0(round(coefs[11,],digits), 
                    ifelse(ps[11,]<0.01, "***", ifelse(ps[11,]<0.05, "**", ifelse(ps[11,]<0.1, "*", "")))),
             paste0("(",round(ses[11,],digits),")"),
             paste0(round(coefs[12,],digits), 
                    ifelse(ps[12,]<0.01, "***", ifelse(ps[12,]<0.05, "**", ifelse(ps[12,]<0.1, "*", "")))),
             paste0("(",round(ses[12,],digits),")"))

xtable(tab)




#####
## Table A9: Robustness to other years
#####

## make a table with three panels

vars <- c('bonw_02_wnind_st', 'bonw_02_wnind.1950_st', 'bonw_02_wnind.1925_st')

all <- as.data.frame(all); fam <- as.data.frame(fam)

coefs <- matrix(NA, 3, 10)
ses <- matrix(NA, 3, 10)
ps <- matrix(NA, 3, 10)
rsq <- matrix(NA, 3, 10)
ns <- matrix(NA, 3, 10)

for (i in 1:length(vars)){
  all$bon_st <- all[,vars[i]]
  fam$bon_st <- fam[,vars[i]]
  # run models
  mod1 <- lm(coup ~ bon_st, all)
  mod2 <- lm(coup ~ bon_st + nind_log, all)
  mod3 <- lm(coup ~ bon_st + biz + pol + mil + syrian + immig + nind_log, all)
  mod3b <- lm(coup ~ bon_st + biz + pol + mil + syrian + immig + nind_log + clust, all)
  mod4 <- lm(coup ~ bon_st, fam)
  mod5 <- lm(coup ~ bon_st + value_log_st + cshare2_st + noconsump, fam)
  mod6 <- lm(coup ~ bon_st + value_log_st + cshare2_st + noconsump +
               ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st, fam)
  mod7 <- lm(coup ~ bon_st + value_log_st + cshare2_st + noconsump +
               ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st +
               biz + pol + mil + syrian + immig, fam)
  mod8 <- lm(coup ~ bon_st + value_log_st + cshare2_st + noconsump +
               ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st +
               pol + mil + syrian + immig + nind_log, fam)
  mod8b <- lm(coup ~ bon_st + value_log_st + cshare2_st + noconsump +
                ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st +
                pol + mil + syrian + immig + nind_log + clust, fam)
  mods <- list(mod4, mod5, mod6, mod7, mod8, mod8b, mod1, mod2, mod3, mod3b)
  for (t in 1:length(mods)){
    mod <- mods[[t]]
    se <- coeftest(mod, vcov = vcovHC, type = "HC1")
    coefs[i,t] <- summary(mod)$coef['bon_st',1]
    ses[i,t] <- se['bon_st',2]
    ps[i,t] <- se['bon_st',4]
    rsq[i,t] <- summary(mod)$r.squared
    ns[i,t] <- length(mod$residuals)
  }
}

digits = 3
tab <- rbind(paste0(round(coefs[1,],digits), 
                    ifelse(ps[1,]<0.01, "***", ifelse(ps[1,]<0.05, "**", ifelse(ps[1,]<0.1, "*", "")))),
             paste0("(",round(ses[1,],digits),")"),
             paste0(round(rsq[1,], digits=digits)),
             paste0(round(ns[1,], digits=1)),
             paste0(round(coefs[2,],digits), 
                    ifelse(ps[2,]<0.01, "***", ifelse(ps[2,]<0.05, "**", ifelse(ps[2,]<0.1, "*", "")))),
             paste0("(",round(ses[2,],digits),")"),
             paste0(round(rsq[2,], digits=digits)),
             paste0(round(ns[2,], digits=1)),
             paste0(round(coefs[3,],digits), 
                    ifelse(ps[3,]<0.01, "***", ifelse(ps[3,]<0.05, "**", ifelse(ps[3,]<0.1, "*", "")))),
             paste0("(",round(ses[3,],digits),")"),
             paste0(round(rsq[3,], digits=digits)),
             paste0(round(ns[3,], digits=1)))
xtable(tab)




#####
## Figure A6: robustness to bonacich measures varying 1/delta w own wealth at 0
#####

## importers data

coef <- paste0('bonw',seq(1,9,1),'_wnind_own0')
for (i in 1:length(coef)){
  t <- fam[,coef[i]]
  t_st <- (t-mean(t, na.rm=T))/sd(t, na.rm=T)
  fam <- cbind.data.frame(fam, t_st)
  setnames(fam, c('t_st'), c(paste0(coef[i], '_st')))    
}

vars <- c(paste0(rev(coef), '_st'))
betas <- matrix(NA, length(vars), 6)

for (i in 1:length(vars)){
  fam$t <- fam[,vars[i]]
  mod <- lm(coup ~ t, data = fam)
  betas[i,1] <- summary(mod)$coef['t',1]
  betas[i,2] <- coeftest(mod, vcov = vcovHC, type = "HC1")['t',2]
}

betas[,3] <- betas[,1] + 1.64*betas[,2]
betas[,4] <- betas[,1] - 1.64*betas[,2]
betas[,5] <- betas[,1] + 1.96*betas[,2]
betas[,6] <- betas[,1] - 1.96*betas[,2]

pdf('03_Figures/coef_centrangeweight_own0.pdf')
par(xpd = F, mfrow = c(1,1), mar = c(5, 5, 4, 4))
n <- length(vars)
plot(x = NULL, y = NULL, xlim = c(0, dim(betas)[1]+1), ylim = c(-.05, .2), xaxt = "n", ylab = 'Coefficient',
     xlab = "", cex.axis = 1.5, cex.lab = 1.5)
points(x = seq(1, dim(betas)[1], 1), y = betas[,1], col = 'firebrick1', pch = 16, cex = 1.3)
for (i in 1:n){
  points(x = rep(i, 2), y = c(betas[i,3], betas[i,4]), 
         col = 'firebrick1', type = "l", lwd = 4)
  points(x = rep(i, 2), y = c(betas[i,5], betas[i,6]), 
         col = 'firebrick1', type = "l", lwd = 2)
}
abline(h = 0, col = 'black', lty = 2)
mtext(expression(paste("Scaling parameter ", frac(1,delta))),
      cex = 1.2, padj = 14)
axis(side = 1, at = seq(1,n,1), padj = .4,
     labels = c(expression(frac(9,10 * lambda)),
                expression(frac(8,10 * lambda)),
                expression(frac(7,10 * lambda)),
                expression(frac(6,10 * lambda)),
                expression(frac(5,10 * lambda)),
                expression(frac(4,10 * lambda)),
                expression(frac(3,10 * lambda)),
                expression(frac(2,10 * lambda)),
                expression(frac(1,10 * lambda))))
dev.off()


## all elite data

coef <- c(paste0('bonw', seq(1, 9, 1), '_wnind_own0'))
for (i in 1:length(coef)){
  t <- all[,coef[i]]
  t_st <- (t-mean(t, na.rm=T))/sd(t, na.rm=T)
  all <- cbind.data.frame(all, t_st)
  setnames(all, c('t_st'), c(paste0(coef[i], '_st')))    
}

vars <- c(paste0(rev(coef), '_st'))
betas <- matrix(NA, length(vars), 6)

for (i in 1:length(vars)){
  all$t <- all[,vars[i]]
  mod <- lm(coup ~ t, data = all)
  betas[i,1] <- summary(mod)$coef['t',1]
  # betas[i,2] <- robust(mod)['t']
  betas[i,2] <- coeftest(mod, vcov = vcovHC, type = "HC1")['t',2]
}

betas[,3] <- betas[,1] + 1.64*betas[,2]
betas[,4] <- betas[,1] - 1.64*betas[,2]
betas[,5] <- betas[,1] + 1.96*betas[,2]
betas[,6] <- betas[,1] - 1.96*betas[,2]

pdf('03_Figures/coef_centrangeweight_own0_all.pdf')
par(xpd = F, mar = c(5, 5, 4, 4))
n <- length(vars)
plot(x = NULL, y = NULL, xlim = c(0, n+1), ylim = c(-.05, .2), xaxt = "n", ylab = 'Coefficient',
     xlab = "", cex.axis = 1.5, cex.lab = 1.5)
points(x = seq(1, n, 1), y = betas[,1], col = 'firebrick1', pch = 16, cex = 1.3)
for (i in 1:n){
  points(x = rep(i, 2), y = c(betas[i,3], betas[i,4]), 
         col = 'firebrick1', type = "l", lwd = 4)
  points(x = rep(i, 2), y = c(betas[i,5], betas[i,6]), 
         col = 'firebrick1', type = "l", lwd = 2)
}
abline(h = 0, col = 'black', lty = 2)
mtext(expression(paste("Scaling parameter ", frac(1,delta))),
      cex = 1.2, padj = 14)
axis(side = 1, at = seq(1,n,1), padj = .4,
     labels = c(expression(frac(9,10 * lambda)),
                expression(frac(8,10 * lambda)),
                expression(frac(7,10 * lambda)),
                expression(frac(6,10 * lambda)),
                expression(frac(5,10 * lambda)),
                expression(frac(4,10 * lambda)),
                expression(frac(3,10 * lambda)),
                expression(frac(2,10 * lambda)),
                expression(frac(1,10 * lambda))))
dev.off()



####
## Figure A5c-d: plot high and low reachability
#####

pdf('03_Figures/hist_reach_all.pdf')
hist(all$reachability, main = "", xlab = "Reachability", cex.lab = 1.8, cex.axis = 1.5)
dev.off()

pdf('03_Figures/hist_reach_imp.pdf')
hist(fam$reachability, main = "", xlab = "Reachability", cex.lab = 1.8, cex.axis = 1.5)
dev.off()




#####
## Table A10: robustness to removing outliers
#####

## recoding 

all$bonw_02_wnind_st_l3 <- car::recode(all$bonw_02_wnind_st, "3:hi=NA")
fam$bonw_02_wnind_st_l3 <- car::recode(fam$bonw_02_wnind_st, "3:hi=NA")
all$bonw_02_wnind_st_l1 <- car::recode(all$bonw_02_wnind_st, "1:hi=NA")
fam$bonw_02_wnind_st_l1 <- car::recode(fam$bonw_02_wnind_st, "1:hi=NA")
all$bonw_02_wnind_st_log <- standardize(log(all$bonw_02_wnind + 1))
fam$bonw_02_wnind_st_log <- standardize(log(fam$bonw_02_wnind + 1))
all$bonw_02_wnind_st_rank <- standardize(rank(all$bonw_02_wnind_st, na.last = 'keep'))
fam$bonw_02_wnind_st_rank <- standardize(rank(fam$bonw_02_wnind_st, na.last = 'keep'))
all$bonw_02_wnind_st_dec <- standardize(as.numeric(all$bonw_02_wnind_dec))
fam$bonw_02_wnind_st_dec <- standardize(as.numeric(fam$bonw_02_wnind_dec))


vars <- c('bonw_02_wnind_st', 'bonw_02_wnind_st_l3',
          'bonw_02_wnind_st_l1', 'bonw_02_wnind_st_log',
          'bonw_02_wnind_st_rank', 'bonw_02_wnind_st_dec')

all <- as.data.frame(all); fam <- as.data.frame(fam)

coefs <- matrix(NA, 6, 10)
ses <- matrix(NA, 6, 10)
ps <- matrix(NA, 6, 10)
ns <- matrix(NA, 6, 10)

for (i in 1:length(vars)){
  all$bon_st <- all[,vars[i]]
  fam$bon_st <- fam[,vars[i]]
  # run models
  mod1 <- lm(coup ~ bon_st, all)
  mod2 <- lm(coup ~ bon_st + nind_log, all)
  mod3 <- lm(coup ~ bon_st + biz + pol + mil + syrian + immig + nind_log, all)
  mod3b <- lm(coup ~ bon_st + biz + pol + mil + syrian + immig + nind_log + clust, all)
  mod4 <- lm(coup ~ bon_st, fam)
  mod5 <- lm(coup ~ bon_st + value_log_st + cshare2_st + noconsump, fam)
  mod6 <- lm(coup ~ bon_st + value_log_st + cshare2_st + noconsump +
               ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st, fam)
  mod7 <- lm(coup ~ bon_st + value_log_st + cshare2_st + noconsump +
               ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st +
               biz + pol + mil + syrian + immig, fam)
  mod8 <- lm(coup ~ bon_st + value_log_st + cshare2_st + noconsump +
               ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st +
               pol + mil + syrian + immig + nind_log, fam)
  mod8b <- lm(coup ~ bon_st + value_log_st + cshare2_st + noconsump +
                ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st +
                pol + mil + syrian + immig + nind_log + clust, fam)
  mods <- list(mod4, mod5, mod6, mod7, mod8, mod8b, mod1, mod2, mod3, mod3b)
  for (t in 1:length(mods)){
    mod <- mods[[t]]
    se <- coeftest(mod, vcov = vcovHC, type = "HC1")
    coefs[i,t] <- summary(mod)$coef['bon_st',1]
    ses[i,t] <- se['bon_st',2]
    ps[i,t] <- se['bon_st',4]
    ns[i,t] <- dim(mods[[t]]$model)[1]
  }
}

digits = 2
tab <- rbind(paste0(round(coefs[1,],digits), 
                    ifelse(ps[1,]<0.01, "***", ifelse(ps[1,]<0.05, "**", ifelse(ps[1,]<0.1, "*", "")))),
             paste0("(",round(ses[1,],digits),")"),
             round(ns[1,],digits),
             paste0(round(coefs[2,],digits), 
                    ifelse(ps[2,]<0.01, "***", ifelse(ps[2,]<0.05, "**", ifelse(ps[2,]<0.1, "*", "")))),
             paste0("(",round(ses[2,],digits),")"),
             round(ns[2,],digits),
             paste0(round(coefs[3,],digits), 
                    ifelse(ps[3,]<0.01, "***", ifelse(ps[3,]<0.05, "**", ifelse(ps[3,]<0.1, "*", "")))),
             paste0("(",round(ses[3,],digits),")"),
             round(ns[3,],digits),
             paste0(round(coefs[4,],digits), 
                    ifelse(ps[4,]<0.01, "***", ifelse(ps[4,]<0.05, "**", ifelse(ps[4,]<0.1, "*", "")))),
             paste0("(",round(ses[4,],digits),")"),
             round(ns[4,],digits),
             paste0(round(coefs[5,],digits), 
                    ifelse(ps[5,]<0.01, "***", ifelse(ps[5,]<0.05, "**", ifelse(ps[5,]<0.1, "*", "")))),
             paste0("(",round(ses[5,],digits),")"),
             round(ns[5,],digits))
xtable(tab)






